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AESTRACT 

Substructuring methods are In common use in structural mechanics problems 
where typically the associated linear systems of algebraic equations are 
positive definite* Here these methods are extended to problems which lead to 
nonpositive definite, nonsymmetric matrices. The extension is based on an 
algorithm which carries out the block Gauss elimination procedure without the 
need for interchanges even when a pivot matrix is singular. Examples are 
provided wherein the method is used in connection with finite element 
solutions of the stationary Stokes equations and the Helmholtz equation, and 
dual methods for second-order elliptic equations. 
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1. THE SUBSTRUCTURING ALGORITHM IN THE POSITIVE DEFINITE CASE 


The use of substructuring techniques in the numerical solution of problems 
governed by positive definite partial differential equations is in widespread 
use. The most notable case is found in structural mechanics, especially in 
connection with the equations of linear elasticity. For the sake of 
simplicity, here we describe the technique for the Dirichlet problem for the 
Poisson equation. Specifically, suppose we want to solve 

-Au = f in ft 

( 1 ) 

u = 0 on 3ft 


where ft is, say, an open bounded region in 1R? with boundary 3ft. We 

subdivide the region ft into open subregions ft^, i = l,*»»,ra, such that 
__ m _ 

ft = U ft. and ft Oft. = 0 for i * j. We denote by r, , 1 < i < j < m 
^ l J ij 

the interfaces between regions and ft ^ , i.e., JTj . Of course, 

for particular choices of i and j in a given subdivision, may be 

empty. A sketch of a particular example with m «* 5 is given in Figure 1. 



Figure 1. A subdivision of a region into five subregions. 



We also subdivide ft into a finite difference or finite element grid 
which in practice is much finer than the above subdivision of fl into m 
subregions. We choose the two subdivisions so that the interfaces 
coincide with edges of the finite difference or finite element cells. The 
discretization of (1) proceeds in the usual manner. The essence of the 
substructuring algorithm is found in the particular choice for the ordering of 
the unknowns and equations, i.e., columns and rows, in the linear system 
resulting from the discretization of (1). Specifically, all unknowns and 
equations associated with the interior of a substructure are numbered 

sequentially, one substructure at a time, and unknowns and equations 
associated with the interfaces are grouped together and numbered last. 

For example, in a typical finite difference discretization of (1), one 
associates equations and unknowns with nodes in the grid. In this case, we 
would group together all the unknowns in subregion together and number 

them first, then proceed to etc., and finally to Then we would 

number all the unknowns along the interfaces T. . 1 _< i, j jC m. The 
equations would be numbered in the same way.* Likewise, in a finite element 
discretization of (1), some unknowns (trial functions) and equations (test 
functions) are associated with nodes or edges and these are 


*The subdivision and numbering method described here applies to difference 
methods with stencils involving only nearest neighbors. The method may be 
extended in an obvious manner, e.g., by defining the interfaces to be more 
than one grid point in thickness, to methods having stencils with a greater 
degree of connectivity. 
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2 

numbered in the same manner as in the finite difference case above. In 
addition there may be test and trial functions more naturally associated with 
the finite elements themselves, and the equations and unknowns associated with 
these functions are grouped together with the other ones associated with the 
interior of the corresponding subregion 

The net result of the above numbering schemes is that the linear system 
resulting from the discretization of (1) has the form 



In (2), the matrices A^, i - l,***,m, in the finite element case, result in 
the case of both the test and trial functions being associated with the 
interior of the subregion i = l,***,m, respectively, while the matrix 
Aq results from test and trial functions associated with the interfaces T ^ , 
1 < i < j < m. The matrices C^, and Bj, represent trial, respectively test, 


2 

Again, the method described here applies to the case where the test and trial 
functions vanish outside the elements which contain the associated node or 
edge. However, by defining the interfaces to be one or more elements thick, 
the method may be easily extended to other cases, e.g., cubic B-spline test 
and trial functions. 



functions associated with the interior of and test, respectively trial, 

functions associated with the interfaces. The vectors U^, i * 1, 
respectively denote the unknowns associated with the interior of 
i = while Ug denotes the unknowns associated with the interfaces. 

All of these associations can also be made in the finite difference case. 

It is well-known that the coefficient matrix of the linear system (2), 

resulting from a discretization of (1), is symmetric and positive definite. 

T T T 

Indeed, - A^, for i = l f *«»,m and Aq - A^. It is also easy to 

see that the matrices Aj, i = l,»*»,m, are themselves positive definite. In 

fact, these matrices are exactly the ones which would result from the 

analogous discretization of the problems 

Au = f in 

(3) 

u= 0 on 3ft^ 

for i - l,***,m, where 3fl^ denotes the boundary of Note that this 

boundary may consist of both interfaces and a portion of the boundary 3Q of 
fl, as is the case for ^4* and ^5 *- n figure 1, or may consist wholly 

of interfaces as is the case for f l ^ * n that figure. Discretization of (3) 
results in a linear system with a coefficient matrix A^, and thus A^ is 
clearly symmetric and positive definite. We note that even in the case of the 
Neumann problem, i.e., the boundary condition in (1) is replaced by 3u/3n = 0 
on 3 ft, the matrices A^ in (2) would still be, at least in the finite 
element case, symmetric and positive definite. This is so because the 
problem (3) associated with the matrix A^ is now given by 
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Au » f in 
■ 0 on 

u ■ 0 on 3fl^Or^, j - l,»»»,ra, 


(4) 


where we have set T^ = T^. Since * 8 never empty, the matrix 
associated with (4) is symmetric and positive definite.^ 

With the matrices A^, i = l,««*,m, being positive definite, one may 
proceed to solve (2) by a block elimination procedure. Symbolically, we may 
express the first m stages of this procedure by the relations 


u i - A I 1(F i - B i V’ i-i, •••,«; (5) 

which uniquely express in terms of data and the interface unknowns Uq. 

The last stage of the process requires the solution of the linear system 


where 


DU 


0 


G 



l C aT 1 B and G = F 
i=l 1 


m 


I 

i=l 




( 6 ) 

(7) 


If on 3^0 3f2 something other than Dirichlet data is specified, then the 
matrix A^ also contains rows and columns associated with test and trial 
functions associated with nodes or edges on that portion of the boundary. 

^Of course, the fact that A^, i = are positive definite may be 
deduced directly from the fact the coefficient matrix of (2) is positive 
definite, i.e., the former is a necessary condition for the latter. 
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Of course, in (5) and (7) the inverses are not explicitly computed, but rather 
appropriate linear systems are solved. The solvability of the system (6) 
follows whenever the system (2) is solvable. In fact, if the system (2) is 
positive definite, so is the matrix D [1]. Once (6) is solved for Uq, (5) 
yields U^, i « l,«**,m. 

Although we have described the substructuring algorithm in the context of 
the Poisson equation, the method can be applied in a similar manner to any 
positive definite problem. As noted above, the method has encountered great 
success in structural mechanics problems. However, in other fields where the 
governing equations are not positive definite or symmetric one may still order 
the equations and unknowns to produce linear systems such as (2), but these 
may not always be solved by a standard block elimination procedure. In the 
next two sections we describe a procedure to solve (2) even in the case of the 
matrices A^ being singular and show how the method may be implemented 
through an elimination procedure. In Section 4 we describe examples which 
lead to singular matrices A^. Finally, in Section 5 we give some concluding 
remarks. 

Incidentally, in almost all situations the use of a properly implemented 

substructuring algorithm will result in savings in computational costs when 

compared to a banded elimination procedure. For example, consider a 

discretization of Poisson^s equation on a unit square. Suppose we have H 

o 

subregions in each direction so that m = M and suppose that each subregion 

is further subdivided by introducing an n x n grid. Thus, there are a total 

4 4 

of Mm points in each direction. Banded elimination requires 0(M n ) 
operations, while the above substructuring algorithm can be implemented in, at 

/ / o 

most, 0(Mn + n n ) operations. We note that this particular problem is not 



particularly well-suited for substructuring methods. Also, the relative 
advantage of substructuring is greater when one considers three-dimensional 
problems or systems of partial differential equations. 

We also note that substructuring ideas in connection with preconditioning 
techniques have been discussed in [2], 


2. THE SOLUTION ALGORITHM IN THE GENERAL CASE 

We begin by describing a method for solving (2) in the case where the 
matrices are singular. The algorithm described here is a special case of 

a more general algorithm which applies to arbitrary matrices with arbitrary 
subdivisions into blocks, e.g., the matrix has no special structure and the 
matrices A^ may not only be singular, but may even be rectangular. The more 
general algorithm is described in [3]. We will describe the algorithm as 
applied to (2) and we will make use of pseudo-inverses in order to simplify 
the initial presentation. However, we emphasize that the algorithm may be 
implemented without the need for the explicit calculation of any pseudo- 
inverses; such an implementation is discussed in the next section. This is 
similar to the observation that the algorithm contained in (5)-(7) may be 
implemented without explicitly computing any inverses, e.g., by solving linear 
systems. 

The system (2) is equivalent to 


A i U i + B i U 0 = V i = !»•••. m. 


c i »i + *0 D o ■ F o- 


( 8 ) 


(9) 



- 8 - 


Now, may be orthogonally decomposed in the form 

U i “ Y i + Z i» 1 = (10) 

where 

Aj Zj « 0, i “ 1 , • • • ,tn, (11) 

and Y^ is orthogonal to all vectors satisfying (11). In particular, 

Yj Z ± = 0, i - 1, • • • ,m. (12) 

Substitution of (10)— (11) into (8) yields that 

A i Y i = F i " B i U 0* 1 = 1 ’’**’ m * (13) 

Since is orthogonal to the null space of (13) yields that 

Y i = A i (F i ” B i V’ i = (14) 

where A* denotes the pseudo-inverse of A i# This relation states that 
is uniquely determined from the data and Uq. Note that (8) yields no 
information concerning Z^ as is to be expected since A^, = 0. 

Substituting (10) and (14) into (9) yields that 

m 

DU q = G - l C i Z ± (15) 


where 



-9- 


D = A„ 


- j. C i B i and G “ F o - j, C i A 1 V 


(16) 


i=l 


i=l 


We may also decompose Uq in the form 


u o - Y o + z o 


(17) 


where 


DZ 0 - 0 


(18) 


and Y 0 is orthogonal to all vectors satisfying (18). In particular, 


*0 Z 0 ' °- 


(19) 


Substitution of (17 )— ( 18) into (15) yields that 

m 

DY = G - l C Z (20) 

i-1 

and, since Yq is orthogonal to the null space of D, (20) yields that 

Y 0 = D + (G - l C Z ). (21) 

i-1 

Again, it is not surprising that (15) yields no information concerning Zq. 
Substitution of (17) and (21) into (14) then yields that 

Y i = A it F i “ B i D+ ^ G " C j Z j)^ " A t B i Z 0 (22) 

for i = 1 , • • • ,m. 
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At this point we have shown that Y^, i = 0,«**,m, may be uniquely 
expressed in terms of Z^, i = 0,»*»,m, by (21) and (22). It remains to show 
how to find the latter. The first step is to multiply (13) by (i - A+). 

Since A^ A* A^ = A^ , we have that 

(I - A £ A^)(F 1 - B i U Q ) =0, i = !*••• f m, 


or substituting (17) and (21), 


(I - A ± A^)[F i - B a Z q - B ± D + (G - ^ C ^ Z.)) =0, i = m. 


Now suppose we are able to determine bases for the null spaces of A^, 
i = l,***,m, and D. We collect each of these basis sets into matrices 
i = (),•*•, m, i.e., N^, i = 0,»**,m, have linearly independent columns, 


DNq = 0 and =0, i = l,***,m, 


(24) 


and the columns of Nq, respectively N^, span the null space of D, 
respectively A^, i = l,«««,m. The number of columns in is, of course, 
the dimension of the corresponding null spaces. Now, we may write that 
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where 


and 


R iO " t 1 - A i A t) B i N 0 - H i " (I - A i A l)[ F i - B i D+ G 1 
R A j = (I - A £ A+)B 1 D + Cj Nj , j = !;•••, in. 


(27) 


Now letting 



H 



(28) 


(26) may be expressed in the form 


RA = H. 


(29) 


In general, R is a rectangular matrix. The number of rows in R is equal to 
the sum of the number of rows of the matrices A^, i = l,***,m, and the number 
of columns of R is equal to the sum of the dimensions of the null spaces 
of A^, i = ].,•••, in; and D. It can be shown [3] that the system (29) is a 
consistent system, and we may find its solution, for example, by forming 

(R T R)A = R T H. (30) 


Suppose we can solve (30) for A. Then (28) yields A^, i = 0,»»»,m, (25) 
then yields Z^, i = 0,**»,m, (21) and (22) yields Y^, i = 0,**»,m, and 
finally (10) and (17) yield the solution U^, i = 0,»»»,m, of (2). 
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The algorithm described here is related in the following manner to the 
block elimination algorithm in Section 1. Suppose that the matrix of (2) and 
all the A^'s and D are nonsingular. Then, the algorithm of this section 
reduces to the standard block Gauss elimination procedure. Indeed, in this 
case, A+ => A^, D + = D * and ■ 0 so that and the latter are 

determined uniquely by (14) and (21). Note the correspondence, in this case, 
between ( 14)— ( 15) and (5)-(6). 

In the more general case, i.e., some or all of the A^'s and D being 

singular, it can be shown [3] that the rank deficiency of (30) is exactly that 

of the original coefficient matrix in (2). Therefore, if the latter is 

nonsingular, then so is R T R and then A in (30) is uniquely determined. 

Since the Z±'s and Y^'s are uniquely determined from A, the algorithm 

produces the unique solution of (2). If the matrix of (2) is singular, so 

is R R and (30) does not have a unique solution. However, (30) may be 

solved anyway, either in terms of arbitrary parameters or by adding 

constraints. The number of parameters or constraints is equal to the 

dimension of the null space of R R which in turn is the same as the 

dimension of the null space of the coefficient matrix in (2). In any case, 

once a particular A is determined, then and are also determined. 

In particular applications to the solution of partial differential 

equations, the dimension of the system (30) is small compared to that of the 

T 

system (2). Indeed, typically dim(R R) = 0(m) , the number of subregions. 

For example, the dimension of the null spaces of the matrices and D may 

T 

be one or zero, in which case dim(R R) < m + 1. 
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3. AN ELIMINATION IMPLEMENTATION 

We begin by restating the algorithm of the previous section. Given the 

matrices , B.,***,B , C.,***,C and the vectors F n> *««,F , we find 

0* nr 1 9 9 nr 1 * 9 m 0 9 m 9 

vectors u o ,## * ,U m satis fy* n g (2) by the following procedure. 

1. Compute A+ and A* B^ for i = l,***,m. 

2. Compute , i = l f *** 9 m 9 whose columns constitute a basis for the null 

space of A^, i = respectively. 

3. Compute C^(a* B^), C^(A* F^) and C^ for i = l,***,m. 

m m 

4. Compute D = A Q - J C^A* B jl ) and G = F Q - J C^A* F i ). 

5. Compute D + G. 

6. Compute Nq whose columns constitute a basis for the null space of D. 

7. Compute D + C^ for i = l,«**,m. 

8. Compute the matrices 

R i0 = B 1 N o " A i ( A ± B i) N o for 1 = 

R ij = V D+ c j N j) " V A i B ±X D+ c j N j) for = i»***. m 

and the vectors 

H ± = ^ - B i (D + G) - A lC Aj F) + AjA* B i )(D + G) for i = l,***,m. 

9. Assemble the results of step 6 into the matrix R and vector H 

according to (28) and then compute R T R and R T H. 

T T 

10. Solve the linear system R RA = R H for A and then compute A^, 
i = 0,***,ra, according to the partition of (28). 

11. Compute A^ for i = 0,***,m. 
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m m 

12. Compute Y - (D G) - ]> C N A = (D G) - £ C Z . . 

i=l 1 i=l 

13. Compute Uq = Y 0 + Z 0* 

14. Compute = (A+ F^ - (A* B^Uq for i = l,**»,m. 

15. Compute = Y i + for * = l,«*»,m. 

Other than steps 1, 2, 5, 6, and 10, the above algorithm requries only matrix 
and matrix-vector multiplications. In this section we show how to carry out 
the other operations required by the algorithm through an elimination 
procedure. In particular, we will not need to explicitly calculate any 
pseudo-inverses of matrices. 

We first describe how to carry out steps 1 and 2. Consider the linear 

system. 

A ± Q = S- (B i ,F i ,0) (31) 

where the right-hand side matrix S consists of the matrix B^, the vector 

F^, and some additional columns of zeroes. The number of these additional 

columns should be greater or equal to the dimension of the null space of 
A^. This dimension will actually be determined during the elimination 

procedure. ^ We now proceed to solve (31) by Gauss elimination with partial 
pivoting. If the matrix A^ is singular, then one or more times during the 
elimination procedure we will not be able to locate a nonzero pivot element. 
In fact, the number of times this occurs is exactly the dimension of the null 
space of A^. However, at such an occurrence, the corresponding column is 

**See Section 5 concerning the effects that roundoff errors may have on the 
determination of this dimension. 



-15- 


already in the eliminated form so that we may skip over to the next column and 
continue the elimination process* At the end of the process, (31) has been 
reduced to the form 

\ Q = J = CV^.O) (32) 

where is upper triangular and in row echelon form. When A^ is 

singular, A^ will have zeros at the pivot location for exactly those columns 
for which no nonvanishing pivot element was found. 

We now proceed to backsolve (32). No difficulty is encountered until a 
row is reached for which the pivot entry of A^ is zero. For the columns 
of Q corresponding to and F^, we may arbitrarily set (to something 

other than zero) the entry in the row corresponding to the zero pivot of A^ . 
Then the backsolve procedure may continue until we reach another zero pivot 
entry, at which time we again arbitrarily specify an entry in the columns of 
Q corresponding to the columns and F^ of S. While all this is going 

on we are also solving (32) for the columns corresponding to the zero columns 
of S. For these columns, whenever a zero pivot entry is encountered in A^ , 
one of the elements in the corresponding row is set to one while the rest are 
set to zero. Each time a zero pivot entry is encountered, a different column 
is chosen for which one sets the arbitrary element to one. At the end of this 
backsolve procedure, (32) yields that 

Q = (L.K,^) . 

A 

Here the columns of form a basis for the null space of A^ and L and 

a 

K are particular solutions of the systems. 
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L = and K = F^,. (33) 

A A 

The final step Is to orthogonalize the columns of L and K with respect 
to the columns of to yield 

rw //v /s/ N 

Q = (L,K,N ± ). 

Since A^N^ = 0, L and K are still solutions of (33). Moreover, the 

columns of L and K are orthogonal to the null space of A^ and, 

therefore, are minimum norm solutions. By the uniqueness of the minimum norm 
solution, we have that 

L = A* B ± and K = A* 

Thus the above elimination procedure has accomplished the tasks of steps 1 and 
2 of the algorithm. 

The tasks of steps 5, 6, and 7 can be accomplished in an analogous 

manner. Also, if the matrix R R is nonsingular, then it may be easily 

solved by an ordinary Gauss elimination procedure. If it is singular then a 
solution in terms of arbitrary parameters may be determined in a manner 
similar to the above procedure for the system (31). We note that any sparsity 
or structure inherent in the matrices may be exploited in the above 

procedure. However, in general, the matrix D will be dense. We will return 
to this point in the concluding section. 
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4. EXAMPLES 

The Stationary Stokes Equation 

Consider the stationary Stokes equations for the slow flow of a viscous 

fluid in a bounded region in IE?. These are given by 

A_u - grad p = _f in fl 

div u = 0 in !1 (34) 

ij “ 0 on 

Here £ denotes the velocity, p the pressure, _f the given body force and 

the viscosity coefficient has been absorbed into p and _f. Clearly, the 

pressure cannot be determined uniquely since we may add an arbitrary constant 

to the pressure and still satisfy (34). 

A finite element approximation of the solution (u,p) of (34) may be 

defined as follows. Given finite-dimensional spaces V* 1 and for the 

ll ll ll 

discrete velocity and pressure fields, we seek €V and p £S such that 
/(grad xx 1 : grad _v^ - p* 1 div v^)dft = -J dft for all v^€V^ 

a n 

(35) 

/ q 1 div u 1 dft = 0 for all 

a 

Here we assume that the elements of satisfy the boundary condition in 

(34), By choosing bases for the spaces V* 1 and S* 1 , (35) can be expressed as 

a linear algebraic system for the coefficients in the basis function 
h l 

expansions of u n and p . 
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u u 

Now it is well-known that arbitrary choices of spaces V n and S may 
not yield stable or accurate solutions. However, there are now known many 
element pairs for which (35) yields optimally accurate solutions [4], [5], 

[6]. One such pair is described as follows. Suppose denotes a 

triangulation of the region ft. We denote by a finer triangulation 

derived from S h by subdividing each triangle in into four congruent 

triangles by joining thp midsides. See Figure 2. We define S* 1 to consist 
of piecewise constant functions over the triangulation 




Figure 2. A triangle in and the corresponding triangles in 

and V* 1 to consist of piecewise linear functions over the triangulation 
which are continuous over ft and vanish on 3ft. This combination is known to 
be stable and be optimally accurate [6].^ The basis functions for V* 1 are 
easily associated with the vertices of the triangulation while the basis 

functions for S* 1 are associated with the triangles in the triangulation S^. 


^See below for the necessary restriction on the pressure which yields this 
result. 
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Now let us consider a substructuring technique for the solution of (35) • 
We assume that the interfaces T ^ between subregions are made up of edges of 
the triangulation so that these interfaces do not cut across pressure 

triangles. One may easily arrange a numbering scheme for the unknowns and 

equations which yields a linear system of the form (2), For example, 
consists of all velocity unknowns associated with vertices of located in 

the interior of the subregion and all pressure unknowns associated with 

the triangles of which are also in Note that Uq contains only 

velocity unknowns, namely those associated with vertices which lie on the 

interfaces but not on 8 fl. 

We have not constrained the pressure space and therefore the system (2) 

corresponding to this discretization of (34) is singular. In fact, its rank 

defficiency is one, and the null vector corresponds to the pressure function 

which is constant over ft. On the other hand, the velocity approximation is 

uniquely determined by (2) [6]. Furthermore, it is easy to see that the 

submatrices A_,***,A are singular. In fact, these matrices are exactly 
1 m 

those which arise from the analogous discretization of the problem. 

Aii - grad p = f_ in ft^ 
div = f in ft^ 
u = 0 on Sft^. 

Thus each of the matrices A^ has a single local pressure null vector, i.e., 
the dimension of N^ is one and N^ corresponds to the pressure function 
which is constant over ft^. On the other hand, since the velocity field can 
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be uniquely determined from (2) and since Uq consists of only velocity 

unknowns, the matrix D in the linear system (15) is nonsingular, i.e., 

Ng b 0. Thus, in this case, the system (30) has dimension m and has a one- 
dimensional null space, the latter following from the fact that the system (2) 
itself has a one-dimensional null space. 

If we choose the pressure space S* 1 to consist of piecewise linear 

functions over the triangulation which are continuous over ft, while 

retaining the same velocity space, the situation changes drastically. For 
example, now the basis functions for S* 1 are more easily associated with the 
vertices of S^. Now contains pressure unknowns corresponding to 

vertices in S, which are in the interior of ft, or lie on 3ft.Pl3ft. More 
n i I 

important, Uq now contains pressure unknowns associated with vertices of 
which lie on T ^ but not on 3ft. In this case the matrices are 

nonsingular and the matrix D is singular with a one-dimensional null space. 


The Helmholtz Equation 

Now consider the problem 

Au + Xu =* f In ft 
u = 0 on 3ft 


(36) 


where X is not near an eigenvalue of the operator -A. Standard finite 
element or finite difference discretizations of (36) yield linear algebraic 
systems with coefficient matrices which are symmetric and indefinite, but 
which certainly may, by using a partial pivoting strategy, be stably 
inverted. Now consider the following * specif ic situation. Let ft be the 
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square (0,ir) x (0,7t) and let X = 13/4, Since the eigenvalues of -A for 

2 2 

this region are given by (n + m ), m,n = 1,2,**«, we see that X = 13/4 is 
not an eigenvalue and therefore the problem (36) leads to nonsingular 
coefficient matrices. Now, suppose we consider solving (36) by using the 

substructuring algorithm with the two subregions 0^ = (0,2 tt/3) x (0,tt) and 
58 (2 tt/3,tt) ft ( 0 , tt ) • Then the matrices in (2) correspond to the 

coefficient matrix for the analogous discretization of the problem 

Au + Xu = f in ft^ 

(37) 

u = 0 on 8ft^. 

2 2 

But the eigenvalues of -A for the region f l ^ are given by (n + 9m /4), 
ra,n = 1,2,3,***, so that X = (13/4) is an eigenvalue of -A for the region 
ft^ and therefore the matrix is singular even though the system (2) is 

not. 

Admittedly, this example is somewhat pathological in the sense that for 
random choices of regions, subregions, and parameters X, the probability is 
zero that the matrices A ^ in (2) will be singular. However, for particular 
choices of X , ft and ft^, one or more of the matrices A^ may be singular; 
after ^11, the above example is not really all that far-fetched. Of course, 
if any of the A^s are singluar, the situation may be remedied by choosing a 
different subdivision of the region 0; this in turn implies a complete 
reassembly of the coefficient matrix in (2). On the other hand, the algorithm 
of Sections 2 and 3 may be used whether or not any of the matrices A^ are 
singluar. 
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There is a small but nonvanishing probability that for some of the 
problems (37) A, although not an eigenvalue of -A for the region ^ , is 
close to such an eigenvalue. If A is close enough to such an eigenvalue, 
the matrix A^, in finite precision arithmetic, may be mistakenly determined 
to be singular by the algorithm of Section 3. However, this will be the case 
only when the difference between A and an eigenvalue is much smaller than 
the discretization error, i.e., of the order of the unit roundoff error of the 
machine, and no serious effect on the accuracy of the solution should result. 


Dual Methods for Second-Order Elliptic Equations 

For a third example, we consider dual methods for second-order elliptic 
partial differential equations. An example of these are methods based on the 
complementary energy principle in linear elasticity. For simplicity, we here 
consider the problem 


and 


u « V<|> in f 2 

div u_ = f in ft 

ii*n_ =0 on r 

<t> = g on r 2 


(38) 


where again = 8ft denotes the boundary of the bounded region ft CUT 

and n denotes the unit outer normal to 8ft. A finite element approximation 
of (38) may be obtained by choosing finite-dimensional spaces V* 1 and S ^ 
and then seeking u^y* 1 an( j such that 
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/(uV + div v^)dn =» f gv h .ndft Vv^gy* 1 
ft r 2 

/ <|> h div u h dn = / f^ h ¥<|> h es h . 

ft ft 

We assume that the elements of V* 1 satisfy the boundary condition on in 

(38). The boundary condition on <f> is natural in this formulation, which is 
one of its advantages. 

In [7], the following choice of V* 1 and S* 1 was shown to yield stable 
and optimally accurate approximations, at least for polygonal domains. First, 
we subdivide ft into quadrilaterals, and then subdivide each quadrilateral 
into four triangles by drawing the diagonals. For V* 1 we take all continuous 
piecewise linear vector fields with respect to the resulting triangulation and 
then define S* 1 = div V* 1 . The resulting space S* 1 can be shown to be a 

subspace of all piecewise constants over the triangulation. See [7] for 
details. 

In the implementation of the substructuring algorithm, we assume that the 
interfaces coincide with some of the edges of the quadrilaterals which 

initially defined our finite element triangulation of ft, i.e., the interfaces 
do not cut through any of these quadrilaterals. The test and trial functions 
from V h are associated with nodes while those from S* 1 are associated with 
the interior of the quadrilaterals. The matrices in (2) now correspond 

to the discretization of the problem 

_u = V<}> and div _u = f in ft^ 
u*n_ =0 on TjOaft^, = g on T^PlSft^ 


( 39 ) 



and 


ji = 0 on 

Because of the last boundary condition, the .problem (39) is over constrained 
insofar as the variable u us concerned. Nevertheless, if « 0, 

i.e., a given subregion does not have part of its boundary coincide with that 
part of 3ft on which data for <J> are given, then the problem (39) can only 
determine <{> to an additive constant. This, for example, would be the case 

for subregion Figure 1, i.e., an interior subregion. For such 

situations, i.e., = 0, the matrix in (2) will again be singular, 

with a one-dimensional null space. Since (38) always uniquely determines u 9 
the matrix D of (16) will be nonsingular. The rank deficiency of the system 
(30) will be one or zero, depending on whether or not Y^ has vanishing 
measure, i.e., whether or not the problem (38) uniquely determines <j>. 


5. CONCLUDING REMARKS 
Determination fo Zero Pivot Elements 

A crucial step in the elimination algorithm presented in Section 3 is the 
determination of when all the elements in a column to be eliminated are 
already zero. This is necessary for the determination of the null spaces of 
the matrices A^ and D. In practice one would declare an element to vanish 
whenever its magnitude is less than some prescribed tolerance which should be 
proportional to the unit roundoff error of the machine. This naturally leaves 
open the possiblity of a very small but nonzero element being mistaken for a 
vanishing element. This situation can be avoided, at least when one is 



-25- 


solving partial differential equations, by first using high enough precision 
arithmetic, e.g., 60 or 64 bit floating point arithmetic, and by making sure 
that the algorithms used are stable. The former is easily arranged, while the 
latter points out the importance of rigorous mathematics. Indeed, if an 
algorithm is stable, as are the ones discussed in Section 4, and the machine 
precision is high enough, one should not encounter nonzero elements which are 
comparable in magnitude to the unit roundoff error unless the matrix in hand 
is singular or very nearly singular. 

An alternative to the use of elimination type procedures is, of course, to 
employ methods based on orthogonal transformations. At the price of greater 
computational expense, such methods are less susceptible to ill effects due to 
roundoff error. 


Parallelism 

One of the attractions of substructuring algorithms is the obvious 
inherent parallelism both in the assembly and solution stages. The sets of 
matrices and vectors * = can each be assembled 
independently. Furthermore, at least in the finite element case, we may write 
the matrix Ag and the vector Fg in the form 


m 

A ° = i=l Aoi> 


m 


l 

i=l 


Oi 


(40) 


where the matrix Ag^ and the vector Fg^ represent the contribution to the 
matrix Ag and vector Fg coming from region Each of the sets (A^, 

F 0i^* * = ma y b e assembled in parallel. Thus, in the assembly stage, 

the sets ^i» B i ,C i ,F i» A 0i ,F 0i^ * * = ma y be assembled In parallel. 
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For example, each of the above sets may be assembled on separate processors, 
with no need for interprocessor communications. At the end of the assembly 
process, the concatenations of (40) must be performed. This step is not 
parallelizable, but represents a minor portion of the assembly process. 

There is also a large degree of parallelism in the solution algorithm 
described at the beginning of Section 3. Steps 1, 2, and 3 are completely 
parallelizable, again with no interprocessor communications necessary. 
Furthermore, if the appropriate information can be transferred to the 
processors, steps 7, 8, 11, 14, and 15 and a portion of step 12 can also be 
computed in parallel. The only relatively major steps which are not 
parallelizable are steps 5 and 6. 

The issue of parallelism in connection with substructuring algorithms has 
been studied in [8] in the context of a specific three-dimensional positive 
definite problem. That paper contains a discussion of operation counts which, 
for the most part, is also relevant in the present context. 

Three-Dimensional Problems 

As pointed out above, the major nonparallel steps in the computation are 
embodied in steps 5 and 6 in the algorithm of Section 3. Even on a serial 
machine these steps may be costly since, in general, they involve dense 
matrices. In two-dimensional problems, by keeping the number of subregions 
relatively small compared to the total number of elements in the 
triangulation, the size of these dense calculations can be kept small, i.e., 
the size of D can be of the order of the square root of the size of the 
A^s. The latter usually are sparse, e.g., banded. A similar arrangement in 
three-dimensional problems would, in general, lead to a matrix D whose size 
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is of the order of the two-thirds power of the size of the A^'s, which may be 
unacceptably large* Furthermore, in steps 1 and 2 of the algorithm, the 
number of right-hand sides would be approximately equal to the number of 
columns of D and the size of the A^"s may be too large, when relatively 
few subregions are used. Therefore, for three-dimensional problems one must 
be especially careful to implement the algorithm in an efficient manner as 
possible. 

These potential difficulties can be mitigated in a variety of ways. For 
example, many of the right-hand sides in the computations of step 1 of the 
algorithm are zero because any column of which corresponds to an 
interface unknown which is not associated with 8ft^ would vanish. The 
corresponding row of is also zero. Thus, one can avoid computations 
involving linear systems with zero right-hand sides and multiplications by 
zero vectors. The savings possible, in storage and computing time, by 
accounting for these features are relatively higher for three-dimensional 
problems. 

Although, in general, the number of interface variables may be large for 
three-dimensional problems, in practice it is often the case that specific 
features of the domain ft lead to a small number of such unknowns. For 
instance, in a wing-fuselage configuration, it is natural to consider the wing 
and fuselage to be different subregions and the interface between these two 
substructures is relatively small in extent. Indeed, it was exactly In this 
type of application that the terminology "substructuring" arose. 

Finally we consider the most serious problem, namely that of the size of 
the matrix D. However, even here a judicious implementation can effect great 
savings. As a simple Illustration consider the subregion structure of Figure 
3 where we have now labeled the interface boundaries by i- !,•••, m-1. 
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Figure 3. An example subdivision of the region ft. 

It is natural to order the interface unknowns Uq one interface at a time, 
e.g., first those on then those on etc. It is not hard to see that 
the matrix D for this example is block tridiagonal, i.e., the unknowns 
corresponding to the interface are connected only to the unknowns on the 
interfaces ^±- 1 , F^, and By taking advantage of features such as 
this, the cost of step 5 and 6 of the algorithm can be greatly reduced, 
especially in three-dimensional settings. We note that these ideas are 
similar to those connected with one-way direction algorithms for positive 
definite problems [9]. 



- 29 - 


REFERENCES 


[1] F. Gantmacher, The Theory of Matrices , Chelsea, New York, 1960. 

[2] G. Golub and D. Mayers, "Use of preconditioning over irregular regions," 
in Computer Methods in Applied Science and Engineering , VI, (R. Glowinski 
and J. Lions, Eds.), 1983, pp. 3-14. 

[3] M. Gunzburger and R. Nicolaides, "Elimination with noninvertible pivots," 
Linear Algebra Appl ., 64, 1985, pp. 183-189. 

[4] V. Girault and P.-A. Raviart, Finite Element Approximation of the Navler- 
Stokes Equations , Springer, Berlin, 1979. 

[5] J. Boland and R. Nicolaides, "Stability of finite elements under 
divergence constraints," SIAM J. Numer. Anal ., 20, 1983, pp. 722-731. 

[6] J. Boland and R. Nicolaides, "Stable and semistable low order finite 
elements for viscous flows," SIAM J. Numer. Anal ., 22, 1985, pp. 474-492. 

[7] G. Fix, M. Gunzburger, and R. Nicolaides, "On mixed finite element 
methods for first-order elliptic systems," Numer. Math., 37, 1981, pp. 


29-48 



- 30 - 


[8] L. Adams and R. Voigt, "A methodology for exploiting 
finite element process,'' in Proceedings of the NATO 
Speed Computations , (J. Kowolik, Ed.), Springer-Verlag, 
373-392. 

[9] A. George and J. Liu, Computer Solution of Large Sparse 


parallelism in the 
Workshop on High 
Berlin, 1984, pp. 

Positive Definite 


Systems , Prentice Hall, Englewood Cliffs, New Jersey, 1981. 



Standard Bibliographic Page 


1. Report No. NASA CR-178122 
ICASE Report No. 86-30 


2. Government Accession No. 


3. Recipient’s Catalog No. 


4. Title and Subtitle 

ON STRUCTURING ALGORITHMS AND SOLUTION TECHNIQUES 
FOR THE NUMERICAL APPROXIMATION OF PARTIAL 
DTFFERF.NTTAT, EQUATIONS 


5. Report Date 

Hay 1986 


6. Performing Organization Code 


7. Author(s) 

M. D. Gunzburger, R. A, Nicolaides 


Performing Organization Report No. 

86-30 


9. Performing Organization Name and Address 

Institute for Computer Applications in Science 
and Engineering 

Mail Stop 132C, NASA Langley Research Center 
Hampton. VA 23665-5225 


10. Work Unit No. 


11 


. Contract or Grant No. 

NASl-18107 


12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, D.C. 20546 


13. Type of Report and Period Covered 
Contractor Report 


14. Sponsoring Agency Code 

- 505-31-83 - 01 


15. Supplementary Notes 

Langley Technical Monitor: 
J. C. South 

Final Report 


Submitted to 
Mathematics 


Applied Numerical 


16. Abstract 


Substructuring methods are in common use in mechanics problems where 
typically the associated linear systems of algebraic equations are positive 
definite. Here these methods ae extended to problems which lead to 
nonpositive definite, nonsymmetric matrices. The extension is based on an 
algorithm which carries out the block Gauss elimination procedure without the 
need for interchanges even when a pivot matrix is singular. Examples are 
provided wherein the method is used in connection with finite element 
solutions of the stationary Stokes equations and the Helmholtz equation, and 
dual methods for second-order elliptic equations. 


17. Key Words (Suggested by Authors(s)) 
substructuring, parallel computing 


18. Distribution Statement 

64 - Numerical Analysis 
Unclassified - unlimited 


19. Security Classif.(of this report) 

20. Security Classif.(of this page) 

21. No. of Pages 

22. Price 

Unclassified 

Unclassified 

32 

A0 3 


For sale by the National Technical Information Service, Springfield, Virginia 22161 
NASA Langley Form 63 (June 1985) 


NASA-Langley, 1986 










End of Document 



